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APPLICATION OF MULTIVARIABLE SEARCH TECHNIQUES 
TO THE OPTIMIZATION OF AIRFOILS IN A 
LOW SPEED NONLINEAR INVISCID FLOW FIELD 

By D. S. Hague and A. W. Merz 
SUMMARY 

Multivariable search techniques are applied to a 
particular class of airfoil optimization problems. These 
are the maximization of lift and the minimization of dis- 
turbance pressure magnitude in an inviscid nonlinear flow 
field. A variety of multivariable search techniques contained 
in an existing non-linear optimization code, AESOP, are 
applied to this design problem. These techniques include 
elementary single parameter perturbation methods, organized 
search such as steepest-descent, quadratic, and Davidon 
methods, randomized procedures, and a generalized search 
acceleration technique. Airfoil design variables are 
seven in number and define perturbations to the profile of 
an existing NACA airfoil. The relative efficiency of the 
techniques are compared. It is shown that elementary one 
palrameter at a time and random techniques compare favorably 
with organized searches in the class of problems considered. 

It is also shown that significant reductions in disturbance 
pressure magnitude can be made while retaining reasonable 
lift coefficient values at low free stream Mach numbers. 
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The optimal solutions reported here were obtained by 
application of a generalized multivariable search code, 

AESOP, orginally constructed under contract to the National 
Aeronautics and Space Administration's Office of Advanced 
Research and Development. Original documentation of this 
code is given in references 1 to 3; an outline of the analysis 
underlying this code is presented below. 


MULTIVARIABLE SEARCH 

The general non-l1near multivariable optimization problem 
is concerned with the maximization or minimization of a pay-off 
or performance function of the form 

(|) = (|)(ct^)ii = l»<4»...»N (1) 

Subject to an array of constraints 

Cj**Cj(cx^)“0» 3®lf2*...fP (2) 
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The are the independent variables whose values are to be 
determined so as to maximize or minimize the performance 
function subject to the constraints of equation (2). 

The O.J may be looked upon as the components of a control 
vector, a. In a space of dimension N. Since maximization 
o:f a function is equivalent to minimization with a change 
of sign, it suffices to discuss the case in which the per- 
.formance function is to be minimized. 


Multivariable optimization problems involving inequality 
<*'?nstraints may also be encountered. If the constraints are 
applied directly to the independent variables 


L H 

a < < a 

i ’ i 


(3) 


the inequality constraints define a region of the control 
space within which the solution must lie. Inequality con- 
straints on functions of the independent variables similarly 
restrict the region in which the optimal solution is to be 
obtained. In this case 

Inequality constraints can be used to restrict the search 
»|*egion directly, or, alternatively, they may be applied in an 
indirect fashion by a transformation to equality constraints. 
Several transformations may be used for this purpose. For 
example, let an equality constraint, C|^, be defined by the 
transformation 

.(Ej- - E, < Ej; 

(5) 

C„ - 0 . eL < e,^ < e« 

'(E« - E^)^ < E^ 



Constraining Cj^ to zero will result in the constraint of 
equation (4) being satisfied. 

Problems involving equality constraints can be treated 
as unconstrained problems by replacing the actual performance 
function* ()>(ai')t by an augmented performance function, 
where 


P 



j = l 


It can be shown that, provided the positive weighting constants 
Uj are sufficiently large in magnitude, minimization of the per 
formance function subject to the constraints , equation (2) , is 
equivalent to minimization of the unconstrained penalized 
performance function defined by equation (6). This approach 
permits search techniques for finding unconstrained minima 
to be applied in the solution of constrained minima problems 
at the cost of some increased complexity in the behavior of 
the performance function, the performance response surface. 

In practical application, the weighting constants U. are 

w 

determined adaptively on the basis of response surface be- 
havior. 

Alternatives to this approach are available, notably 
Bryson's approach to the steepes t-descent search, reference 
4. This method has been exploited in connection with the 
numerical solution of variational problems encountered in 
the optimization of aerospace vehicle flight paths, refer- 
ences 5, 6, and 7. However, the use of such techniques 
implies smoothness of the response surface. This smoothness 
may not be present in general; hence, the AESOP code is 
limited to the less restrictive penalty function approach of 
equation (6) 
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Numerical Solution of Non>L1near Multivariable 
Optimization Problems 

This section is devoted to a discussion of the search 
algorithms for solution of non-linear multivariable optimi- 
zation problems available in the AESOP code. A wide variety 
of search algorithms have been devised for the solution of 
multivariable optimization problems. Many of these algorithms 
are restricted to the solution of linear or quadratic problems. 
Algorithms of this type must be supplemented by more general 
search procedures if generality of solution is sought; for 
engineering problems tend to lead to non-linear formulation 
with the possibility of discontinuities in both the performance 
function response surface and its derivative. Most of the 
searches which prove effective in these problems combine a 
direction generating algorithm, such as steepest-descent , with 
with a one-dimensional search. Distance traversed through 
the control space in the selected direction is measured by a 
step-size, or perturbation parameter, DP. The object of the 
one-dimensional search is to determine the value of DP which 
minimizes the performance function along the chosen ray and to 
establish the corresponding control vector. 

In practice, the diverse nature of non-linear multivariable 
optimization problems leads to the conclusion that no one 
search algorithm can be uniquely described as being the "best" 
in all the situations which may be encountered. Rather, a 
combination of searches, some of which may be of quite elemen- 
tary nature, provides the most reliable and economical conver- 
gence to the optimal solution. 
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One-dimensionali search . Multivariable search problems 
are reduced to one-dimensional problems whenever a search 
algorithm is used to establish a one-to-one correspondence 
between the control vector and a single scalar perturbation 
parameter, (DP). In such a situation 

= a^(DP), i * 1, 2, , , N (7) 

so that equation (1) becomes 

(|> = (J)(a^) = 4>(DP) (8) 

Similarly, the right hand sides of equations (2) and (6) 
become functions of the scalar perturbation parameter. 

The relationship, equation (7), specifies a ray through 
the control space. As noted above, the objective of the 
one-dimensional search along this ray is to locate the value 
of DP which provides the minimum performance function value. 

Numerical search for the one-dimensional minima can be 
carried out in a local fashion, by the Newton-Raphson method, 
for example, or by a global search of the ray throughout the 
feasible region. The localized polynomial approximation is 
appropriate to the terminal convergence phase in a problem 
solution when some knowledge of the extremal's position has 
been accumulated by the preceding portion of the search and the 
problem involves a smooth function. The global search can be 

I 

used to advantage in the opening moves of a search. In the 
early phase of a search the object is tc isolate the approx- 
imlate neighborhood of the minimum performance function 
value as rapidly as possible, usually with little or no 
foreknowledge of the performance function behavior. One 
measure of the effectiveness of a search algorithm in such 
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a situation is the number of evaluations required to locate 
the minimum point to some prespecified accuracy. It can be 
shown that the most effective method of locating the minimum 
point of a general unimodal function is a Fibonaoai search^ 
reference ..4. In this method, the accuracy to which the mini- 
mum is to be located along the perturbation parameter axis 
must be selected prior to the commencement of the search. 

Since the accuracy required is highly dependent on the behavior 
of the performance function, this quantity is difficult to 
prespecify. 

Prespecification of the accuracy to which the extremal's 
position is to be located can be avoided for little loss in 
search efficiency by use of an alternative search based on the 
so-called golden section, references. This is the method 
employed in the AESOP code one-dimensional search procedure. 
Search by the golden section commences with the evaluation 
of the performance function at each end of the search interval 
and at 6 = 2/(1 + \/5) = .6180339887 of the interval from both of 
these bounding points. This is illustrated in Figure 1. 

The boundary point furthest from the lowest resulting 
performance function value is discarded. The three re- 
maining points are retained, and the search continues in a 
region which is diminished in size by G. The internal point 
at which the performance function is known in the reduced 
interval will be at a distance G of the reduced interval 
from the remaining bounding point of the original interval 
for (1 - G) = 6^. The search can, therefore, be continued 
ini the reduced interval with a single additional evaluation 
of the performance function. It follows after Q evaluations 
of the performance function that the position of the extremal 
point will be known to within R of the original search region 
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Berfomanee Ftmetlon Value 



where 


R » (9) 

To reduce the interval of uncertainty to .00001 of 
the original search interval, about 27 evaluations of the 
performance function are required. For a reasonable number 
of evaluations of the performance function this type of 
search is almost as efficient as a Fibonacci search. 

It should be noted that search by the golden section 
proceeds under the assumption of unimodality; hence, it will 
often fail to detect the presence of more than one minimum 
when the performance function is multimodal. If more than 
one minimum does exist, the one located depends on perfor- 
mance function behavior within the original search interval. 

Multiple Extremals on One-Dimensional Ra y. The one- 
dimensional section search described above is unable to 
distinguish one local extremal from another; it will merely 
find one local extremal. This difficulty can be largely 
eliminated by the addition of some logic to the search, at 
least for moderately well behaved performance functions; 
that is, for functions having a limited number of extremals 
in the control space region of interest. An effective method 
for detecting multiple extremals is to combine the one- 
dimensional search with a random one-dimensional search on 
the same ray through the control space. This is illustrated 
in figures 2 and 3. In figure 2 the response contours of 
a performance function having two minima are illustrated , 
together with the initial points used in a global one- I 
dimensional search by the golden section method. The 
behavior of the function at these points is shown in 
figure 3, The left hand minimum is not apparent from 
these points. If a single random point is added in the 
Interval Lq. the probability of this point revealing the 
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Perforoanee Function Value 



presence of the second mlnlmuin is 

P, ■ L,/L^ (10) 

for any point in the interval AB indicates the presence of 
a local minimum somewhere in the interval AB, and any point 
in the interval BC indicates the presence of a local maximum 
somewhere in the interval BC. In this latter case, there 
must be a minimum of the function both to the left and to the 
right of the newly introduced point. 

If R random uniformly distributed points are added in the 
interval Lot the probability of locating the presence of 
the second minimum becomes 

Pr = 1.0 - (1.0 - 1.,/Lj,)'' (11) 

The function (L|/Lq) is a measure of the performance 
function behavior. For a given value of this behavior function 
the number of random points which must be added to the one- 
dimensional search to provide a given probability of locating 
a second minimum can be determined. 

j The presence' of multiple minima on a one-dimensional cut 
through an N-dimensional space does not necessarily indicate 
that the performance function possesses more than one minimum 
in| a..mul ti-dimensional sense. It may be that the performance 
function is merely non-convex. This is illustrated by figure 
4.1 The performance function behavior on the one-dimensional 
search in figures 2 and 4 is identical. In figure 2 this 
Indicates the presence of two local extremals; in figure 4, 
a non-convex performance function. 

When a one-dimensional search detects the presence of 
multiple extremals in the local sense above, a decision must 
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be made as to which of the apparent extremals is to be 
pursued during the remainder of the search. Here, without 
foreknowledge of the performance function behavior, logic 
must suffice. Typically, the left or right hand extremal, 
the extremal which results in the best performance, or 

even a random choice may be made. 

0 

It should be noted that logic of this type is not cur- 
rently available in the AESOP code. The AESOP one-dimensional 
search procedure has three distinctive phases. First, each 
search algorithm defines an initial perturbation using either 
pajst perturbation stepsize information or a perturbation mag- 
niltude prediction as in the quadratic search below. Second, 
a perturbation stepsize doubling procedure is employed until 
a point exhibiting diminishing performance is generated. 

Third, having coarsely defined the one-dimensional extremal 
position from steps one and/or two, a golden section search 
is employed to locate the extremal with reasonable precision. 


Multiple extremals - general procedure . The multiple 
extremal search technique included in AESOP is based on 
tdpologioally invariant warping of the performance response 
surface. The response surface is warped in a manner which 
retains all the surface extremals but alters their relative 
locations and regions of influence. The region of influence 
of an extremal is defined as the hull or collection of all 
points which lead to the extremal if a gradient path is 
followed. Reducing the region of influence of an extremal 
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decreases the probability of locating a point in the neigh- 
borhood of the extremal if points are chosen at random. 
Again, in an organized multivariable search, the probability 
0 ^ locating an extremal having a small region of ,^i nf 1 uence 
is less than that of locating an extremal having a large 
region of influence. For example, suppose the extremals 
olf the one-dimensional function of figure 5 are to be deter- 
mined in the range < a < by the sectioning approach. 
The four initial values employed in this technique are 
denoted by f^ to f^. 


Following evaluation at these four points, f 4 is dis- 
carded, and the function is evaluated at f 5 » Atthis point 
the right-hand extremal, e2» has been eliminated from the 
search which now inevitably proceeds to the left hand extre- 
mal at ei . 


To find the second extremal, the Function F is warped 
by writing 


i ■ (ciH-a*) [— 2— 2L"| + ct*j a > a’* 


r * T 2 N 

(a*-aL) [ - ■ I + a*; a* > o 

[a*-aL J 


( 13 ) 


where N is a positive integer, and a* is the Vocation of the 
left hand extremal . 

i A typical relationship between K and a is shown in 
figure 6 for the case N = 1. Differentiation of equation 
(13) with respect to a when N = 1 results in 


2[a-a*] . 

-.T ^ n. . j 

Ca^-a*] 


a ^ a’ 


C = ; a < a 

[a*-aL] 


(14) 


Note that as a -► a*, 0 from both the left and 

right. At a = otL and at a = a^, C' = 2. In the regions 
< a < a* and a* < a < a^, i varies parabol ical ly 
with a. Figure 6 illustrates these points. It can be 
seen that a region Aa-| centered about a* transforms into 
a smaller region AC] located in the neighborhood of 
C a*. On the other hand, a region Ao 2 situated in the 
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neighborhood of the upper search limit, maps into a 
wider region in the neighborhood of C = an* If' general, 
the slopes at a = and a = are given by 2N; the 
greater N’, the greater the warping becomes. 


The effect of introducing a moderate warping trans- 
formation on the function of figure 5 is shown in figure 
7. It can be seen from figure 7 that the region of influ 
ence of e^ is reduced, and the region of influence of 62 
increased. On the warped surface search by sectioning 
commences with evaluations of performance at f^ to f^. 
Following these initial evaluations f| is discarded (as 
opposed to the discard of f 4 on the unwarped surface), 
and the function is evaluated at the additional point fg. 
The points f^ and fg straddle the extremal e 2 which is 
now inevitably located by further sectioning evaluations. 
Figures 8a and 8b illustrate the warping transfo'rmations 
for N = 1 and N = 10 when the transformation is applied 
at the point a* = .5, the symmetric case. It 
can be seen that when N = 1 , twenty per cent of the 
warped control space corresponds to approximately 45 per 
cent of the unwarped control space in the vicinity of 
the transformation origin (a = .5). When N = 10 twenty 
per cent of the warped control space transforms into 
ninety per cent of the unwarped control space. 


Sectioning Parallel to the Axes . The independent 
variable perturbation algorithm in the sectioning search 
is 

* 0 , i f* r 

« DP, i = r r « 1 , 2, . . . , N (15) 
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TRANSFORMED FUNCTION WITH TWO 
EXTREMALS 
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FIGURES 8(a), to 8(d). WARPING TRANSFORMATION N « 1 to 4 
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FIGURES 8(i) and 8(j). WARPING TRANSFORMATION N * 9 and 10 
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This is simply the parametric or univariate search 
approach. All but one of the independent variables are held 
cc vstant while a one-dimensional search parallel to the 
variable axis determines the best value of the remaining 
variable, a^. The variable is then set to this value, 
and the process is repeated with one of the remaining 
independent variables. When all N independent variables 
have been perturbed in this way, a sectioning search cycle 
has been completed. 

The N-dimensi onal search can then be continued with 
another cycle of sectioning or by one of the other search 
techniques described below. In practice, it has been found 
advantageous to perturb the independent variables in a 
random order Within each sectioning cycle. The method can 
be used in conjunction with either a local or a global 
search as outlined in the two preceding sections. The 
behavior of this search in the solution of a straightforward 
two-variable optimization problem is illustrated in Figure 9 
It may be noted that the AESOP code searches from boundary 
to boundary in each variable using a golden section search 
procedure. 

Sectioning to Define Local Sensitivities . The sec- 
tioning search can readily be applied to the problem of 
performance or constraint sensitivity determination. Thus, 
by the device of omitting the updating of each control 
variable following the sectioning search on the r^*^ 
parameter, the sequence of sectioning searches is performed 
about a fixed nominal point. When such a search is per- 
formed in the vicinity of a known extremal point, the 
penalties for off-optimal design can be assessed. Away 
from an extremal point, the search merely provides local 
sensitivities in a similar manner to the manual perturbation 
methods employed in conventional trial and error design 
evolution. 


23 





i 

\ 


I 

1 




Steepes t-Descent Search . The steepest-descent search 
algorithm is 

(Aa) - [|£]^K3]-’{K2)| 

■S /US DC rK3T‘{DC} 

Ki- K2 [K3>*tK2l 

- [W]-> C||]^ [Kj]-* {DO (16) 

Here, the matrix W is the metric tensor of the control space 


and serves to 

def i ne 

a generalized 

measure for the magnitude 

of a control 

vector 

perturbation. 

The vectors {3(j>/3o} and 

{3C/3a> are defined 

as 




3<f) 


3^1 




3ai 

* 3012* 

3a 1 
n 



and 1 

3C 

3C 

3C 



1 

da-| 

• 3ag • 

* * 3a 

n 




respectively. The K matrices are defined as 


''i = LIf J cw] ■■ (|^ } 

(17) 

{K2}= [|^][W]-{|±} 

(18) 

[K3> [|§][W3- [|£)T 

(19) 
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The perturbation parameter, (DP), is defined by 

top) = LAaJ [W3{io) (.20) 

The vector is desired change in the constraint 
functions. For an unconstrained problem, (16) reduces to 

(ao) • -CW]^{|i) (21) 

The performance function change associated with the pertur- 
bation of equation (16) is 

0* - - U^2j[K3]'*{K2)y^^(DP)^- LDC| CK3]" (DC) 

+ Ll^2j[K3]'‘ {DC) (22) 

Equation (16) does not specify a one-dimensional 

! 

search directly since the perturbation parameter (DP) 
and each component of the constraint vector change DC 
can be independently specified. This difficulty is 
conveniently eliminated if the components of DC are 
expressed in terms of the perturbation parameter. Let 
(DP) and DC be arbi trari ly assigned, say (DPq) and D^q, 
respectively. Now consider the one parameter set of 
values for DC defined by 

^ = (§^)- ^0 ( 23 ) 

It follows from equations (16) and (22) that (23) spec- 
ifies a one parameter family of perturbations in which 
the non-linear performance and constraint functions vary 
linearly with (DP), to the first order. 


Equations (16) to (22) are valid for small pertur- 
bations in the independent variables provided the der- 
ivatives involved are continuous in the region of the 
control space defined by equation (20). In practice, 
When this condition is not satisfied, the steepest- 
djescent algorithm can be used to locate a promising 
direction for a one-dimensional search provided the 
cjerivatives are computed numerically. In this case, 
however, equation (22) ceases to provide an accurate 
indication of performance function behavior along the 
specified ray. 


When dealing with performance and constraint 
functions having continuous first derivatives, the 
perturbation parameter value to be used in equation (16) 
can be determined from a second order Taylor’s expan- 
sion of the performance function behavior in terms of 
DP. The coefficients in this series expansion can be 
readily obtained from the conditions of zero change for 
DP = 0, linear slope for DP = 0, and from the actual value 
of the performance function at a point in the neighborhood 
of the point at DP = 0. This method for determining the 
best perturbation parameter value is discussed in some 
detail in references 5 and 6. When dealing with less reg- 
ular functions, the one-dimensional search by sectioning 
can be used to determine the perturbation parameter value. 
This is the technique employed in the optimization program, 
AESOP, references 1 and 2; for the AESOP code converts all 
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cpnstrained optimization problems to unconstrained problems 
by the penalty function device, equation (6). The resulting 
response surface combines both performance function and 
weighted constraint functions. Inevitably, this surface has 
a more complex topology than that of the unconstrained per- 
formance function. Program AESOP is also limited to the 
penalty function approach to constrained optimization, and, 
hence, it utilizes the reduced algorithm of equation (21) 
rather than the explicit constraint algorithm of equation (16) 


Steepest-Descent Weighting Matrices . The weighting 
matrix introduced in equations (16) and (20) must be positive 
definite to assure a positive distance between any two non- 
coincident points in the control space. Apart from this 
restriction, the choice of weighting matrix is arbitrary. 
Inspection of equation (16) reveals that any descending 
direction is a steepest-descent path for some choice of 
the weighting matrix W. This can be simply illustrated 
when only two independent variables are involved. Figure 10 
depicts a small region of the control space . The per- 
formance function response contours appear as a series of 
parallel lines on this microscopic region of the control 
space. The perturbation zones corresponding to three 
weighting matrix choices are shown. The first zone corres- 
ponds to the choice of a unit matrix for W. It follows from 
equation (20) that for a given value of (DP)^ the search 
zone is a circle of radius (DP). The steepest-descent 
direction is that in which the performance improvement is 
greatest. This is the direction of a line from the origin 

V 
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of the circular search zone to that point on its circumference 
which provides the smallest value of the performance function 
With this choice of weighting matrix, the steepest- 
descent direction is perpendicular to the response contours. 
Paths of this type are illustrated in figure 11 by the solid 
lines emanating from points A and B. From the nominal point 
A, search perpendicular to the performance response contours 
is very efficient. From point B, however, this type of 
search results in the meandering path illustrated. It is 
assumed here that once a steepest-descent direction is located, 
an exhaustive search for the minimum in that direction will 


bei undertaken in view of the high cost of recomputing the 
derivatives in many problems. Even if this were not the 
case, search normal to the response contours can often be 
improved upon. For example, it is obvious that even in the 
straightforward two-dimensional problem of figure 11 the 
dashed search direction is superior. This direction requires 
a priori knowledge of the extremal's position, information 
not normally available. 


Returning to figure 10, the second search zone 
depicted corresponds to the choice of a diagonal matrix 
fdr W. The positive-definite constraint on W requires 
that all diagonal elements of the weighting matrix be 
pcisitive. In this case the search zone becomes elliptical 
with the major and minor axes of the ellipse being parallel 

t 

to the coordinate axes. It may be noted that as either of 
the diagonal elements of W becomes large in relation to the 
remaining element, the corresponding element in W inverse 
together with the predicted change in the associated inde- 
pendent variable becomes small. In the limit this reduces 
the search to a one-dimensional search in the remaining 
coordinate. The perturbation zone then becomes a slit 
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nCURB ll.^STEKrEST-DESCSTJr SEARCH 


parallel to that coordinate axis of length 2 *. (DP), as 
Illustrated in figure TO. In the case il lustrated » the 
steepest-descent path is in the descending direction. 

Finally, the search zone corresponding to the choice 
of an arbitrary positive-definite weighting matrix is 
shown. From equation (20) and the positive-definite con- 
straint on W, the search zone remains elliptical, but the 
principle axes may now have an arbitrary orientation to 
the axes of and a 2 . It follows that since the elliptic 
search zone can have any orientation and eccentricity, any 
direction in the control space is a possible steepest- 
descent path; for in all cases, the path of steepest-descent 
lies in the direction of a line joining the search zone 
origin to the lower point of tangency between the boundary 
of the search zone and the performance function response 
contours. The discussion above may readily be extended 
to control spaces of higher dimensionality. 

When attempting the solution of optimization problems 
by the steepest-descent method, the analyst is constantly 
faced with the problem of choosing a satisfactory weighting 
matrix for the search continuation. The problem is com- 
pounded by the fact that the slopes of the performance 
« 

function with respect to the independent variables can, 
and frequently do, vary by many orders of magnitude. The 
arbitrary choice of a unit matrix in such situations can 
lead to distressingly slow convergence of the numerical 
search; for it is in the nature of many problems that in 
those directions in which the slopes are greatest the re- 
sponse surface is highly non-linear , Only small pertur- 
bations will be successful in the direction of these strong 
control variables. In those directions in which the slopes 
are small, the contours are often relatively linear, and 
targe perturbations may be required in these weak control 
variables. 


In such situations the local steepest-descent direction 
for (W) = (I) is misleading. With this choice of weighting 
matrix perturbations are in proportion to the response 
surface partial derivatives. However, the best direction 
in which to proceed may involve large perturbations in the 
weak control variables of small slope. This behavior is 
illustrated for a two-dimensional case in figure 11 by the 
dashed line emanating from B. 


The problem of choosing a satisfactory weighting matrix 
also arises when the steepest-descent search is applied in 
its variational form, reference 5, and when a combination of 
continuous control variables and parameters are encountered 
a$ in the optimization of multiple-arc problems in flight 
path optimization problems, reference 6. In these references 
it is suggested that the weighting matrices be based on the 
first derivatives of the unconstrained performance function 
with respect to the control. This approach can be used in 
the solution of multivariable optimization problems also, by 
writing 

“ij ° ^ • < = j 

w 

= 0, i j 


In practice, alternate use of the resulting combined weighting 
matrix and the unit matrix tends to provide a reasonable con- 
vergence rate at points well removed from the extremal. The 
AESOP code employs such a matrix in combination with a search 
range non-dimensional i zation term and a learning factor. The 
learning factor emphasizes perturbation of control parameters 
which change in a monotonic direction and de-emphasizes those 
perturbations which fluctuate in sign. 

Random Ray Search . The difficulty of defining a 
suitable control variable metric tensor together with 
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the fact that any descending path Is a steepest-descent 
direction for some choice of metric tensor suggests the 
possibility of searching along a random ray through the 
control space. The algorithm for random ray search Is 

« R^(±0P), 1 « 1, 2, . . .,N (25) 

where the .proportional to the direction cosines of the 
ray, are uniformly distributed random numbers satisfying 

-1.0 < R^ < +1.0, 1 = 1, 2, . . ., N 

The positive sign in equation (25) is taken If is 

negative; the negative sign Is taken when this derivative 
is positive. The method is equivalent to a steepest- 
descent search using a randomly generated metric tensor. 


Quadratic search . An alternative systematic approach 
to, the definition of an arbitrary or empirical weighting 
matrix is provided by second order or quadratic method. It 
can be shown, for example, in reference 1, that on an 
elliptic second order response surface the weighting 
matrix 




(26) 


vAien xised in the steepest-desoent method will iitinediately define 
the optimal point 


{a*> ■ (aQ) + {M 


( 27 ) 


^ Also known as the Newton-Raphson method 
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where {6a} is computed from equation (21) with (DP)* = .5K^. 
On a more general non-linear response surface, equation (27) 
merely defines a direction for subsequent search in the 
manner of the steepest-descent technique. This is illus- 
trated in figure 12. Here, the approximating elliptical 
contours computed at point 0 define an approximate extremal 
location at P through equations (26) and (27). Subsequent 
search along the ray OP results in the definition of a one- 
dimensional extremal. This point is then used to fit another 
approximating elliptic contour, and the process is repeated 
until the extremal point at Q is located. 

The quadratic search procedure can be quite rapid in 
control spaces of low dimensionality. In high order spaces 
the approach is usually impractical as a result of the 
requirement to establish the second order weighting matrix 
of equation (26). In many practical engineering problems 
these derivatives cannot be obtained in closed form; in such 
cases the derivatives must be obtained numerically, for 
example, reference 1. Computation of these derivatives 
requires at least (N+l)(N+2)/2 evaluations of 4> at each 
point where an approximating quadratic is employed. Clearly, 
for large N this computation may become impractical in 
computational time. 

Davidon or FI etcher-Powel 1 Method . Davidon's method 
is a hybrid first order/second order technique. The objec- 
tive of Davidon's method is to arrive at a reasonable 
approximation to the second order weighting matrix of 
equation (26) without the use of (N+l)(N+2)/2 evaluations 
of It can be shown that on a quadratic (second order) 
response surface N steepest-descent searches performed in 
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FIG. 12b . ■» QUADRATIC SSARCH ON A HIGSSR ORDER SUBFACH 
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the manner described previously will lead to definition of 
the weighting matrix of equation (26), if the following 
formula Is employed: 



[w],^, - [W], + [A], + [B], 

(28) 

where 

■ [AL- 

(29) 



(30) 


[W]f^ = [I] 

(31) 


Here, L^aJ-| Is the change In position during the 1^^ one- 
dimensional search and 


L* ■ ISJ, 

is the change in gradient vector between the beginning and 
end of the i^ one-dimensional search. On a numerically 
well-behaved function this technique may work well. It 
will find the optimum of an elliptic quadratic function in 
N successive searches. When appreciable numerical noise 
is present in the calculation, or when the one-dimensional 
extremal along the ray is hot defined with precision, the 
method may produce erratic convergence, or convergence 
failure. 

Pattern Search . In the present report, pattern search 
refers to a search which exploits a gross direction revealed 
by one of the other searches. The search algorithm is 

Aa^ ” ^ ® N (32) 

37 



where and aj are the components of the control vector 
before and after the use of a preceding search technique. 

This type of search is illustrated in figure 13 following 
a section search. The combination of a section search 
and a pattern search in the problem illustrated leads 
directly to the neighborhood of the extremal. Repeated 
sectioning, on the other hand, would be a very slowly 
converging process due to the orientation of the contours 
with respect to the axes of the independent variables. 

It may be noted that a simple rotation of the independent 
variable axes by 45® results in sectioning alone becoming 
a: rapidly converging process in this example. The pattern 
search can also be used to accelerate the steepest-descent 
process provided it follows two successive descents as in 
figure 14. 

Adaptive Creeping Search . Adaptive creeping search is 
ai form of small scale sectioning; however, instead of locating 
the position of the one-dimensional extremal on each section 
parallel to a coordinate axis, the coordiate is merely perturbed 

by a small amount, Aa^, in the descending direction. 

The search commences with a small perturbation in one 
Ojf the independent variables', a^; a positive perturbation 
is first made; if this fails to produce a performance 

improvement, then a negative perturbation is tried. If 
neither of the perturbations produces an improved perfor- 
mance value, the variable retains its nominal value, and j 

1 . 

Aa^ is halved. If a favorable perturbation is found, the 
Variable a^, is set to this value, and Aa^ is doubled. The 
process is repeated for each independent variable in turn* 
the order in which the variables are perturbed being 
chosen randomly. At this point an adaptive search cycle 
is complete, and the cycle is then repeated. A two- 
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dimensional illustration of this search is presented in 
figure 15. In the particular problem illustrated, the 
method converges rapidly reaching the neighborhood of the 
extremal within six evaluations. 

The search algorithm can be written in the form 

. (pp) (33) 

whiere S^. is the number of cycles in which the search has 
sujccessf ul 1 y perturbed the r^*^ independent variable, and 
is the number of cycles in which a perturbation of the 
rt:h variable has proved unsuccessful. While this search 
can be looked upon as a one-dimensional approach, this 
viewpoint is somewhat artificial. Here, the scalar quantity 
(DP) merely defines an initial perturbation for each inde- 
pdndent variable. Once started the search proceeds inevi- 
tably to its conclusion, the perturbation in each independent 
variable being adaptively determined according to equation 
(33) on the basis of the performance function response 
contour behavior encountered during the particular problem 
solution. This search can be quite efficient when used in 
combination with the pattern search acceleration procedure. 

Magnification . When studying discrete models of con- 
tinuous systems of the type encountered in certain engineering 
problems such as aerodynamic shaping or structural design 
problems, there is a tendency on the part of some search 
algorithms to achieve a favorable shape before satisfying 
the desired constraint levels. In such cases, when it is 
known that the unconstrained extremal is the null vector, 
a simple magnification search can lead to rapid convergence 

to the desired solution. The magnification algorithm is 

• 

Aa^ = * (DP)» 1=1,2, . . ., N (34) 
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Here (DP) is positive and all components of the control 
vector are to be simultaneously perturbed. Generally, 
the unconstrained extremal point corresponds to the null 
vector; this method may prove efficient. 

Arbitrary Ray Search . In practical design optimization 
a search along an arbitrary multidimensional ray can be of 
utility. For example, when two minimal extv&mal solutions 
appear to be possible, a search on the ray connecting the 
two points should reveal the presence of a maximal extremal 
somewhere on the ray between the two minimal extremals. The 

algorithm for this search is 

2 1 

Aa^ = (a^. - a^.) (DP), i = 1, 2, . . ., N (35) 

where aj- and are the two minimal extremal points. In 
general, and may be any two points in the control space 

Random Point Search . A straightforward Monte-Carlo 
search which examines point designs distributed in a uniform 
random manner within the feasible region is often of utility 
when the response surface is of a complex nature. Such a 
search is included in the AESOP code primarily for use as a 
nominal point design generation procedure. 
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AIRFOIL OPTIMIZATION RESULTS 


Results have been obtained for a number of different low 
speed (M = .1) aerodynamic shaping problems of practical 
interest. Potential flow analysis is based on solution of the 
two-dimensional potential flow equation 

(a^ - u^) 0 + (a^ - v^) 0 - 2uv 0 =0 

' ' '^xx ' ' '^yy '^xy 

where 0 is the velocity potential, u and v are the velocity 

components 

u = 0 , V = 0 

and a is the local speed of sound determined from the energy 
equation and the stagnation speed of sound 


= a^^ - (u^ + v^) 

Solutions are obtained by Jameson's finite difference scheme, 
reference 10. 


The study used all of the numerical search methods described 
above. The nominal airfoil configuration was the NACA 64-206 
airfoil. This airfoil was subject to modifications which en- 
hance its aerodynamic capabilities. The modification was the 
sum of two components: 


a. A continuous binomial additional thickness distribution, 
applied to the upper surface, of the form, 

6Y^(x) = A X ^(1 - X) ^ (O^x^ 1) 

Where A, and are variable parameters to be 
optimized by multivariable search. The quantity A 
is given by 


A 


y(^l ^2) ^ 
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where y is the ma^cimum additional thickness to the 
airfoil. 


b. An additional camber distribution of the form, 

.) 

I \ 

5y^(x) 

(x^< x< 1) 


be selected by multivariable search. These parameters 
denote the location and geometric form of the camber 
of both upper and lower surfaces. 


Xcll - 



Xctl -1 

f * - 
. 1 - 

and 

are \ 


The seven independent parameters of these equations can 
be used to generate a variety of curves, many of which, however, 
are impractical for airfoils. Nevertheless, it is important 
to note that the thickness and camber distributions are both 
’’smooth”, so that the slope is continuous everywhere except 
at the leading and trailing edges of the airfoil. Also, it 
will be found that certain of the parameters can be varied 
only through rather small numerical ranges, in order to retain 
important practical features of the airfoil. 

Results are presented below in three sections, corresponding 
to the following specific applications: 

1. Unconstrained optimization of high-camber airfoil, 

2. Comparative search optimization of low-camber airfoil, and 

3. Constrained optimization of low-camber airfoil. 

The principal purpose of the present study is a demonstra- 
tion of the versatility and relative efficiency of various 
multivariable search options of the AESOP code in the two- 
dimensional airfoil shaping problem. The applications can 
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therefore be regarded as representative, rather than exhaustive. 
This also means that the results obtained can be easily extended 
to include other optimization criteria or other values of the 
reference airfoil or shaping parameters., 

High-Camber Unconstrained Airfoil Optimization 

To demonstrate the AESOP code in the optimization of the 
general airfoil shaping problem, a seven-parameter modification 
to the NACA 64-206 airfoil was tested. The performance cri- 
terion was the lift coefficient of the airfoil, and all of the 
shaping parameters were allowed to vary for this purpose. 

Each of the parameters was given an initial value and 
extreme upper and lower values, as listed in Table I below. 



Parameter 

Minimum 

Initial 

Maximum 


7 

0 

.06 

.09* 

Thickness 





(upper) 


.25 

1.0 

1.0* 


' ^2 

1.0* 

1.0 

3.0 


||n||||| 

.1 

.25 

.9* 

Camber 


0 

.005 

.04* 

(upper 6 





lower) 


1.0 

2.0 

3.0* 


KlH 

1.0 

2.0 

3.0* 


TABLE I. INPUT PARAMETER VALUES IN LIFT MAXIMIZATION 
(♦Signifies optimal parameter value.) 

! 

The optimization procedure was limited to 50 iterations using 
uniform random ray and pattern searches. No other constraints 
were placed on the aerodynamic or geometric properties of the 
airfoil. 
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The results of the optimization study can be presented 
both graphically and numerically. The final, optimal airfoil 
and its associated pressure distribution are shown in Figure 16, 
and the terminal values of the parameters are all at their 
maximum values, except for 82 ' which is at its minimum allowable 
value (see Table I) . 

The following observations can be made concerning these 
results: 

1. The thickness and camber parameters have been combined 
to give a configuration resembling a flapped airfoil, 
which is a standard method of increasing the lift 
coefficient. 

2. The shape of the airfoil is such that separation would 
occur near the trailing edge of the upper surface if 
leading edge stall did not occur first. The pressure and 
lift coefficients are therefore somewhat unrealistic, and 
would presumably require refinement using a viscous theory. 

3. The methods of optimization used are the random ray 
and pattern search methods, and no improvements in 
lift were obtained after the 19th iteration, for all 
parameters are then at their bounding values. The 
initial lift coefficient of 1.347 is increased to 
2.770, as shown in Figure 17(a). This is the limiting 
theoretical lift coefficient for the parameter range 
studied. 

4. l|he adverse pitching moment increases with the lift. 

As shown in Figure 17(b), this variation is nearly 
linear, as was found in earlier optimization studies 
(References 7 through 9). 

For this example, it is seen that practically all performance 
improvements occur between iterations 11 and 18. This type of 
response will be characteristic of unconstrained optimization 
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FIGURE 17(a). LIFT MAXIMIZATION, NO CONSTRAINTS 
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problems which are solved using effective numerical methods. 
Another conclusion of this study is that the constraints can 
comprise ari important part of the input data. Thus, in this 
example, the peak pressure coefficient changed from -5.84 to 
”11.01, but if it had been constrained to a more realistic 
level (greater than -4,0) a very different airfoil contour 
would have resulted. 

Low-Camber Comparative Search 


This portion of the study served to compare the different 
methods of optimizing airfoil performance. For brevity the 
comparison was done with the emphasis on the upper-surface 
parameters, y, and Two nominal configurations were 

examined, corresponding to approximately 3% and 6% additional 
thickness variations. The nominal values of the seven para- 
meters, and the corresponding calculated aerodynamic character 
istics of the resulting airfoil at M = .1 and a = 6°, are 
given in Table II. The geometric appearance of the airfoil, 

and the pressure distribution about the airfoil are shown 
for the 6% modification in Figure 18 using the parameter 


values y = .06, ~ 1-0. 

' The following representative optimization criteria vjere 
chosen for the present brief study: 


a. Minimize peak pressure 


max 


b. Maximize lift (Cj^) 


c. Maximize lift (Cj^) for a given moment (C^) 


Additional input data relates to the optimization 
methods to be used, the number of iterations to be computed, 
and the tolerances permitted on any constraints. 
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Parameter 

Mod. 1 

Mod. 2 

y 

.06 

.03 


i.O 

.75 


1*0 

1.0 


.25 

.25 


.005 

.005 

^3 

2.0, 

1.75 

: £4 

2.0 

2.0 


1.347 

1.086 

! "m 

-.168 

-.094 

*^max 

-5.844 

-5.191 


TABLE II. NOMINAL AIRFOIL CHARACTERISTICS 
a. Minimize Peak Pressure 

Minimization of the high pressure peak at the leading 
edge of the reference airfoil at M = .1 and a = 6°, is a 
good illustrative performance criterion, because small 
geometric airfoil changes can cause large variations in this 
pressure value. Results were obtained using 11 different 
optimization methods, all of which cycled through 100 
iterations on airfoil modification 1 of Table II. The 
results of these calculations are presented in Figure 19, 
which show the improvements in peak pressure as a function 
of iteration number. Also shown in these Figures are the 
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final airfoil shape and pressure distributions. It is seen 
that several methods converge quickly to the minimum "optimal" 
value, while other methods perform poorly. The eleven methods 
are ranked in order of final pressure coefficient in Table III, 
which sums up the final values of Figure 19. Notice that the 
best results follow from randomized and one parameter at a 
time search methods, while "steepest descent" and "quadratic" 
methods give very poor results. It should be noted that the 
computational times include compiling several subroutines, 
program reassembly, loading the combined Jameson/ASEOP codes 
and computer generated plotting. Computation times would be 
smaller if an absolute program element were employed, plots 
omitted and only the relevant parts of AESOP employed. Never- 
theless, the computational times quoted serve to measure the 
relative effectiveness of the eleven searches. 


Method 

m 


Iterations 

Comp. Units 
for 100 
Iterations 

Uniform Random Ray 


32 

1105 

Creeper 


64 

1029 

*Davidon 

1.77 

98 

1083 

Directed Random Ray 

1.85 

24 

1156 

Sectioning (12 Evaluations 
on ID Search) 

1.86 

52 

1351 

*Steepest Descent (Variable W) 

1.96 

81 

1112 

Sectioning (6 Evaluations 
on ID Search) 

2.00 

42 

966 

Monte Carlo 


13 

290^^^ 

♦Jacobson (Homogeneous Functions) 

4.09 

87 

682 

♦Steepest Descent (W = I) 

5.40 

98 

875 

♦♦Quadratic 

5.74 

59 

728 

♦ Methods using 1st derivatives 


• 


♦♦ Method using 1st and 2nd derivatives 



♦** Only ran 15 airfoils 



- 


TABLE III. RELATIVE PERFORMANCE OF OPTIMIZATION METHODS- 


53 







20 N 40 


60 


80 


1 



FIGU?£ 19(d). directed RANDOM RAY METHOD 













The large variations in these "optimal" airfoil shapes 
emphasize that certain numerical methods are much better than 
others in this type of application, and that speed of conver- 
gence should be a major consideration in the choice of numerical 
algorithm. The ranking of the results is roughly in the same 
order as the computational cost, as measured in the last 
column of Table III. Computation times can be significantly 
reduced by elimination of the coding associated with searches 
not used in a given calculation. This was not done in the 
present study however. It can also be noted (Reference 9) that 
the five best methods developed airfoils which have two equal 
pressure peaks on the upper surface, but which are otherwise 
noticeably different. In cases where two peaks are developed the 
forward sharp pressure peak tends to be a result of the potential 
flow analysis employed. In real flows this forward peak is modified 
by viscous effects. /Once an airfoil having two equal pressure peaks 
is developed by the searches further progress in the Cp minimization 
becomes more difficult for only those airfoil perturbations which 
simultaneously reduce both peaks provide further improvement. 

b. Maximize Lift 

The lift coefficient was again chosen as a performance 
criterion for airfoil modification 2, and four optimization 
methods were compared in a 50 iteration test. The thickness 
parameter y was allowed to vary over the interval 0 ^ y ^ .06 , 
and the exponents for the leading and trailing edges were 
required to satisfy .25-^e^:^ 1.0, and 1.0^e2i3.0, respectively* 

The absence of camber parameters excludes the possibility 
of developing a "flap-like" airfoil having a very high lift 
coefficient. 

The results of the optimization methods are shown in 
Figure 20, and it is seen that the methods produce essentially 
equal lift coefficients. All of the methods (directed random 
ray, uniform random ray, steepest descent, and creeper) converged 
to the maximum additional thickness after 20 to 30 iterations, 
after which no further gains in lift were generated for limiting 
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FIGURE 20(a). DIRECTED RANDOM RAY METHOD 



FIGURE 20 (b) . UNIFROM RANDOM RAY METHOD 
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FIGURE 20(C). CREEPER METHOD 



10 N 20 30 40 


I 


FIGURE 20(d). STEEPEST-DESCENT METHOD 
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values of the exponents were attained. These values were y = .06, 
= ^2 - 1-0. The final values of lift# moment and pressure 
coefficients are as given in Table IV. 


Method 

c. 

<=1. 

C 

m 


Cp 

max 


Directed Random Ray 

1.347 

-.168 

5.842 

Uniform Random Ray 

1.342 

-.165 

5.572 

Steepest Descent 

1.347 

-.168 

5.843 

Creeper 

1.347 

-.168 

5.843 


TABLE IV. LIFT MAXIMIZATION BY FOUR NUMERICAL METHODS 

c. Maximize Lift for a Given Moment 

The incorporation of a constraint function into the 
optimization problem adds realism to the problem statement. 

In the present case, the desired moment coefficient was chosen 
at a moderate but representative value, = -.1, and 
airfoil modification 1 was taken as the nominal airfoil. The 
significant initial airfoil parameters and associated lift, 
moment and pressure coefficients have been given in Table II. 

The additional thickness parameters were then adjusted 
using a combination of pattern and random ray techniques, 
with results as summarized in Figure 21. The actual compu- 
tation was carried out for the unconstrained maximization of 
the performance index, 

SO that the degree of violation of the constraint can be 
controlled by the choice of the scalar, k. In the present 
case, it was given the initial value 10,000. This constraint 
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penalty weighting function value is then internally adjusted 
in the AESOP code. 

The final values of the significant parameters were: 

y = .0867 
= .895 

62 = 1-678 

all of which are in the mid-range of the allowable values of 
the parameters. As shown in Figure 21, the moment constraint 
initially causes the lift to decrease, after which a modest 
steady increase in lift is achieved. The initial airfoil 
shape and pressure distributions have been shown in Figure 18, 
and the final results are given in Figure 22, which shows 
the strong sensitivity of the pressure distribution to small 
changes in the parameters. The final values of the lift, 
moment and peak pressure coefficients in this example are: 

= 1.255 

Xj 

C„ = --105 



and here again as in some previous solutions it is noted 
that the two upper surface pressure peaks are equal. 
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FIGURE 22* 64-206 MAXIMIZE LIFT COEFFICIENT 

WITH CM = -.1 



CONCLUSION 


The ability of multivariable search procedures in the 
solution of two-dimensional non-linear inviscid low speed 
aerodynamic shaping problems has been examined. In general the 
non-linear inviscid optimization problem yields results as 
readily as the linear inviscid problem. Reference 11. Elementary 
one parameter at a time and random techniques appear to achieve 
better results than organized searches which require derivatives. 
Only one organized search out of those employed comes near to 
the elementary searches in solution of the present problems. 

This organized search is the first-order Davidon method. 

It is clear from these results that unconstrained lift 
maximization is primarily accomplished by modifications to the 
airfoil trailing edge. In particular this is accomplished by 
the generation of "flap-like” airfoils. 

When pressure peak magnitude is reduced the more successful 
searches define an airfoil which has two equal pressure peaks. 

The first peak lies near the leading edge. The second peak is 
30% to 40% aft of the leading edge and is of a broader more 
gradual nature than that at the leading edge. The forward peak 
is both sharper and narrower, .and would be considerably 
attenuated by the effects of viscosity. Subject to the limitations 
of the viscous theory, however, the pressure optimization process 
yields a logical conclusion. 

It is clear that significant modifications can be made to 
an existing airfoil using a few carefully selected shaping para- 
meters. At low speeds three parameters appear to be sufficient 
for upper surface shaping, seven parameters suffice for both 
upper and lower surface shaping. Use of a small number of 
shaping parameters for two dimensional airfoil shaping indicates 
that successful three dimensional shaping is quite feasible at 
the present time. 
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